Exploring a ZDR Offset ValueΒΆ

import xarray as xr
import pyart
import numpy as np
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import warnings
import xcollection as xc
import glob
import matplotlib.pyplot as plt
import pandas as pd
import cartopy.crs as ccrs
from distributed import Client, LocalCluster
warnings.filterwarnings("ignore")
hv.extension('bokeh')

Read in the dataΒΆ

def create_dataset_from_sweep(radar, sweep=0, field=None):
    """
    Creates an xarray.Dataset from sweeps
    """
    
    # Grab the range
    range = radar.range['data']
    
    # Grab the elevation
    elevation = radar.get_elevation(sweep)
    
    # Grab the azimuth
    azimuth = radar.get_azimuth(sweep)
    
    # Grab the x, y, z values
    x, y, z = radar.get_gate_x_y_z(sweep)
    
    # Grab the lat, lon, and elevation
    lat, lon, alt = radar.get_gate_lat_lon_alt(sweep)
    
    # Add the fields
    field_dict = {}
    for field in list(radar.fields):
        field_dict[field]=(["azimuth", "range"], radar.get_field(sweep, field))
    
    ds = xr.Dataset(
        data_vars=field_dict,
         coords=dict(
             x=(["azimuth", "range"], x),
             y=(["azimuth", "range"], y),
             z=(["azimuth", "range"], z),
             lat=(["azimuth", "range"], lat),
             lon=(["azimuth", "range"], lon),
             alt=(["azimuth", "range"], alt),
             azimuth=azimuth,
             elevation=elevation,
             sweep=np.array([sweep]),
             range=range),
    )
    
    return ds.chunk()
def convert_to_xradar(radar):
    """
    Converts from radar to xradar
    """
    
    ds_list = []
    sweeps = radar.sweep_number['data']
    for sweep in sweeps:
        ds_list.append(create_dataset_from_sweep(radar, sweep))
    
    # Convert the numpy array to a list of strings
    dict_keys = [x for x in list(sweeps.astype(str))]
    
    # Zip the keys and datasets and turn into a dictionary
    dict_of_dsets = dict(zip(dict_keys, ds_list))
    
    return xc.Collection(dict_of_dsets)

List files availableΒΆ

We grab the files from the following directory, after ordering from the data portal and unzipping the archive

radar_files = sorted(glob.glob('../../data/x-band/ppi/*6_PPI.nc'))

Spin up a Dask ClusterΒΆ

cluster = LocalCluster()
client = Client(cluster)
client
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-p297zaac', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-w0kfqa11', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-t_qx4ojw', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-wpkibo54', purging

Client

Client-3d51261e-abc0-11ec-b6c4-acde48001122

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:51851/status

Cluster Info

Loop through and read in the dataΒΆ

ds = convert_to_xradar(pyart.io.read(radar_files[0]))['0']

Clean the DataΒΆ

We will drop a couple of coordinates, and make sure that we don’t include the fill values (which are large negative numbers)

ds = ds.drop(['elevation', 'sweep'])

Set which fields to plot in the interactive plot

fields = ['DBZ', 'VEL', 'ZDR']
for field in fields:
    da = ds[field]
    ds[field] = ds[field].where(da > -999)
ds.compute().hvplot.quadmesh(x='lon',
                             y='lat',
                             ylabel='Latitude',
                             xlabel='Longitude',
                             z=['DBZ', 'VEL', 'ZDR'],
                             rasterize=True,
                             width=600,
                             height=400,
                             geo=True,
                             alpha=0.1,
                             selection_color={'NaN': (0, 0, 0, 0)},
                             #projection=ccrs.PlateCarree(),
                             cmap='pyart_Carbone42',
                             features={'rivers':'110m'},
                            tiles='StamenTerrainRetina').opts("Image", alpha=0.8)

Investigate HistogramsΒΆ

ReflectivityΒΆ

dbz_hist = ds.DBZ.hvplot.hist(bins=np.arange(-40, 40, 1), title='DBZ Values')
dbz_hist

NCPΒΆ

ncp_hist = ds.NCP.hvplot.hist(bins=np.arange(0, 1.01, 0.01), title='NCP Values')
ncp_hist

ZDRΒΆ

zdr_hist = ds.ZDR.hvplot.hist(bins=np.arange(-3, 1.1, .1), title='ZDR Values')
zdr_hist

Plot a 2D histogram using matplotlibΒΆ

ds_subset = ds.where((ds.RHOHV > 0.95) & (ds.NCP > 0.5))
bins = np.linspace(-15, 25, 71)
zdr_bins = np.linspace(-1, 0, 71)
hist = np.histogram2d(ds_subset.ZDR.values.flatten(), ds_subset.DBZ.values.flatten(), bins=70, range=((-1, 0), (-15, 25)))[0]
hist = np.ma.masked_where(hist == 0, hist)

# Create a 1-1 comparison
#x, y = np.meshgrid((-5, 1),(-6, 20))
c = plt.pcolormesh(bins,
                   zdr_bins,
                   hist,
                   cmap='pyart_HomeyerRainbow')

# Add a colorbar and labels
plt.colorbar(c, label='$log_{10}$ counts')
plt.xlabel('DBZ')
plt.ylabel('ZDR')

plt.show()
../_images/zdr-offset-exploration_25_0.png

Plot a 2D Hexbin VisualizationΒΆ

zdr = ds_subset.ZDR.values.flatten()
dbz = ds_subset.DBZ.values.flatten()
df = pd.DataFrame({'ZDR':zdr,
                   'DBZ':dbz})
min_zdr = -2
max_zdr = 1
df_range = df[(df.ZDR >= min_zdr) & 
              (df.ZDR < max_zdr)]
df_range.hvplot.hexbin('DBZ',
                 'ZDR',
                 rasterize=True,
                 cmap='pyart_HomeyerRainbow',
                 logz=True,
                 )

Plot the Filtered DatasetΒΆ

def plot_field(field='DBZ',
               color_bounds=(-20, 40),
               cmap='pyart_Carbone42',
               ds_to_plot=ds_subset):
    plot = ds_to_plot[field].compute().hvplot.quadmesh(x='lon',
                                                        y='lat',
                                                        ylabel='Latitude',
                                                        xlabel='Longitude',
                                                        clim=color_bounds,
                                                        rasterize=True,
                                                        xlim=(-107.1, -106.8),
                                                        ylim=(38.8, 39.1),
                                                        width=600,
                                                        height=400,
                                                        geo=True,
                                                        alpha=0.1,
                                                        projection=ccrs.PlateCarree(),
                                                        cmap=cmap,
                                                        features={'rivers':'110m'},
                                                        tiles='StamenTerrainRetina').opts("Image",
                                                                                          alpha=0.8)
    
    return plot
ref_plot = plot_field('DBZ')
vel_plot = plot_field('VEL', (-15, 15), cmap='pyart_balance')
zdr_plot = plot_field('ZDR', (-2, 1))
(ref_plot + vel_plot + zdr_plot).cols(1)